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The problem of DNA—DNA interaction mediated by divalent counterions is studied using a gener¬ 
alized Grand-canonical Monte-Carlo simulation for a system of two salts. The effect of the divalent 
counterion size on the condensation behavior of the DNA bundle is investigated. Experimentally, 
it is known that multivalent counterions have strong effect on the DNA condensation phenomenon. 

While tri- and tetra-valent counterions are shown to easily condense free DNA molecules in solution 
into toroidal bundles, the situation with divalent counterions are not as clear cut. Some divalent 
counterions like Mg^^ are not able to condense free DNA molecules in solution, while some like 
Mn+2 can condense them into disorder bundles. In restricted environment such as in two dimen¬ 
sional system or inside viral capsid, Mg"*"^ can have strong effect and able to condense them, but the 
condensation varies qualitatively with different system, different coions. It has been suggested that 
divalent counterions can induce attraction between DNA molecules but the strength of the attrac¬ 
tion is not strong enough to condense free DNA in solution. However, if the configuration entropy of 
DNA is restricted, these attractions are enough to cause appreciable effects. The variations among 
different divalent salts might be due to the hydration effect of the divalent counterions. In this 
paper, we try to understand this variation using a very simple parameter, the size of the divalent 
counterions. We investigate how divalent counterions with different sizes can leads to varying qual¬ 
itative behavior of DNA condensation in restricted environments. Additionally a Grand canonical 
Monte-Carlo method for simulation of systems with two different salts is presented in detail. 

PACS numbers: 87.14.gk,87.19.xb,87.16.A- 


I. INTRODUCTION 

The problem of DNA condensation in the presence of 
multivalent counterions has seen a strong revival of inter¬ 
est in recent years. This is because of the need to develop 
effective ways of gene delivery for the rapidly growing 
field of genetic therapy. DNA viruses such as bacterio¬ 
phages provide excellent study candidates for this pur¬ 
pose. One can package genomic DNA into viruses, then 
deliver and release the molecule into targeted individual 
cells. Recently there is a large biophysics literature ded¬ 
icated to the problem of DNA condensation (packaging 
and ejection) inside bacteriophages (for a review, see Ref. 

s. 

Because DNA is a strongly charged molecule in aque¬ 
ous solution, electrostatics and the screening condition of 
the solution play an important role in the structure and 
functions of DNA systems. Specifically, the condensation 
of DNA molecules is strongly influenced by the counte¬ 
rion valence [H- While tri- and tetra-valent counte¬ 
rions are shown to easily condense free DNA molecules 
in solution into toroidal bundles, the situation with di¬ 
valent counterions are not as clear cut. Some divalent 
counterions like Mg+^ are not able to condense free DNA 
molecules in solution, while some like Mn'*'^ can condense 
them into disorder bundles. Similarly, strong electro¬ 
static effect is also observed for DNA condensation in a 


restricted environment such as inside a viral capsid. By 
varying the salinity of solution, one can vary the amount 
of DNA ejected from viruses. Interestingly, monova¬ 
lent counterions such as Na'*'^ have negligible effect on 
the DNA ejection process In contrast, multivalent 
counterions (Z—ions for short) such as Mg+^, CoHex+^, 
Spb"*"^ or Spm'^'^ exert strong and non-monotonic ef¬ 
fects 0. There is an optimal counterion concentration, 
cz.o, where the least DNA genome is ejected from the 
phages. For counterion concentration, cz, higher or lower 
than this optimal concentration, more DNA is ejected 
from phages. The case of divalent counterions is more 
marginal. The non-monotonicity is observed for MgS 04 
salt but not for MgCl 2 salt up to the concentration of 
lOOmM. Such ion specificity for the case of divalent salts 
also present in condensation of DNA in free solution. Qj. 

The non-monotonic influence of multivalent counteri¬ 
ons on DNA ejection from viruses is expected to have 
the same physical origin as the phenomenon of reen¬ 
trant DNA condensation in free solution in the presence 
of counterions of tri-, tetra- and higher valence [Ml- 
Although, divalent counterions are known to condense 
DNA only partially in free solution ii , DNA virus 
provides a unique experimental setup. The constraint of 
the viral capsid strongly eliminates configurational en- 
tropic cost of packaging DNA. This allows divalent coun¬ 
terions to influence DNA condensation similar to that 
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of trivalent/tetravalent counterions. Indeed, DNA con¬ 
densation by divalent counterions has also been observed 
in another environment where DNA configuration is con¬ 
strained, namely the condensation of DNA in two dimen¬ 
sional systems 13|. For virus systems, theoretical fitting 
suggests that the DNA is neutralized at cz,o ~ 75mM for 
divalent counterions, and the short-range DNA attrac¬ 
tion at this concentration is —0.004 A:bT per nucleotide 
base BUI- 

In this paper, we study the problem of DNA conden¬ 
sation in the presence of divalent counterions using com¬ 
puter simulations. The simulation method developed by 
our groups in Ref. B and B is used, expanded and the 
influence of the ion size on the strength of DNA— DNA 
interaction mediated by divalent counterions is investi¬ 
gated The Grand Canonical Monte Carlo simula¬ 
tion for a system of two salts is presented in detail. The 
electrostatic contribution to the free energy of packaging 
DNA into bundles is calculated from simulation. It is 
shown that, if only the non-specific electrostatic contri¬ 
bution is included, divalent counterions can indeed in¬ 
duce DNA reentrant condensation like those observed 
for higher counterion valences. However, correlations 
among divalent counterions are not strong enough to de- 
condense DNA bundles. As already mentioned, experi¬ 
mental results also show that there is a ion specific effect. 
As a first step taken to study this ion specific effect, the 
DNA—DNA effective interaction is calculated from simu¬ 
lation for three different counterion sizes. It is shown that 
varying counterion sizes can have significant impact on 
DNA condensation pictures, which can explained some 
variations among DNA condensation experiments with 
Mg2-|-, or Mn2-|- counterions. 

The paper is organized as follows. In Sec. im the 
Grand-canonical Monte-Carlo is formulated to simulate 
a system of two salts (a divalent salts and a fixed mono¬ 
valent salt from buffer solution). In Sec. Illli the model 
of our system and various physical parameters used in 
the simulation are presented in details. In Sec. ITVl the 
results are presented and their relevance to available ex¬ 
perimental data is discussed. We conclude in Sec. 


II. GRAND CANONICAL MONTE-CARLO 
SIMULATION FOR MIXTURE OF TWO SALTS 


In practical situation, the DNA bundle is in equilib¬ 
rium with a water solution containing free mobile ions 
at given concentrations. Therefore we simulate the sys¬ 
tem using Grand Canonical Monte-Carlo (GCMC) sim¬ 
ulation. The number of ions is not constant during the 
simulation. Instead their chemical potentials are fixed. 
These chemical potentials are chosen in advance by sim¬ 
ulating a DNA—free salt solution and adjusting them so 
that the solution has the correct ion concentrations. An¬ 
other factor that complicates the simulation of DNA con¬ 
densation phenomenon arises from the fact that there are 
both monovalent and divalent salts in solution in experi¬ 


ments. At very low concentration of divalent counterions, 
czi DNA is screened mostly by monovalent counterions. 
To properly simulate the DNA bundle at this low cz 
limit, and to properly capture the screening of electro¬ 
static interactions among divalent counterions by mono¬ 
valent ones, both salts are included in the simulations 


To simulate two different salts present in our system, 
the standard GCMC method for ionic solution B is gen¬ 
eralized to simulate of a system containing a mixture of 
both multivalent and monovalent salts. For simplicity, 
we assume both salts have the same coion (for example, 

Cl“). Thus, a state i of the system is characterized by 
the locations of Niz multivalent counterions, W_|_ mono¬ 
valent counterions and coions. In the grand canon¬ 
ical ensemble of unlabeled particles, the probability of 
such state is given by 

TTi = 4 ■3Ar,+ .3Ar,_ + H+N,+ + - /SUi] 

Az ' A_^ A_ 

(1) 

Here, Z is the grand canonical partition function, /3 = 
l/ksT, ^-z,+,- = h/y'2TTmz,+ -kBT, Ui is the interac¬ 
tion energy of the state i, and fj,z,+,- are the chemical 
potentials of the multivalent counterions, of the monova¬ 
lent counterions and of the coions respectively. 

In a Monte Carlo simulation, a Markov chain of sys¬ 
tem states i is generated with a limiting probability dis¬ 
tribution proportional to tt^. This chain is defined by 
a probability py of transitions from state i to state j. 

A sufficient condition for the Markov chain to have the 
correct limiting distribution is: 


Pij _ ^ 
Pji '^i 


( 2 ) 


As usual, at each step of the chain, a “trial” move to 
change the system from state i to state j is attempted 
with probability qij and is accepted with probability fij. 
Clearly, 


Pij Qij fV 


( 3 ) 


It is convenient to regard the simulation box as consisting 
of V discrete sites {V is very large). Then for a trial move 
where lya particles of species a are added to the system: 


Qij — 


1 

V^. 


( 4 ) 


Conversely, if Va particles of species a are removed from 
the system: 


Qij — 


{Ng - I^a)\ 

NgW- 


( 5 ) 


Putting everything together, equations 0-® give us 
a recipe to calculate the Metropolis acceptance probabil¬ 
ity of a particle insertion/deletion move in GCMC sim¬ 
ulation. For example, if in a transition from state i to 
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state j, a multivalent salt molecule (one Z—ion and Z 
coions) is added to the system, the Metropolis probabil¬ 
ity of acceptance of such move can be chosen as: 


Vice versa, for a “trial” move where one Z—ion is re¬ 
moved from the system and Z monovalent counterions 
are added to the system. 


/m = min{l, jijifji) (6) 


where 


/d 

fji 


__ 

{J^iZ + + Z') 


exp[/3(C/j 


with 


(7) 


yZ+i 

Bz = eXp(/3^2 , (S) 


and 


^z,salt — ^^z + Z^_ (9) 


is the combined chemical potential of a multivalent salt 
molecule. 

On the other hand, if a multivalent salt molecule (one 
Z—ion and Z coions) is removed from the system. 


fij 

fji 


N,zN,--{N,- 

Bz 


Z+l) 


exp[/3(t7i 


( 10 ) 


Similarly, for addition a monovalent salt molecule (one 
monovalent counterion and one coion) in transition from 
state i to state j. 


fij 

fji 


_ BzNjz _ 

Bi {Ni-\. + l)...(Vi_|_ -|- Z) 


exp[,d(C/i 


[/,)]. (16) 


Note that because the system maintains charge neu¬ 
trality in all particle addition/deletion moves, instead of 
using 3 different chemical potentials, to simu¬ 

late the system, only two combined chemical potentials, 
l^zsalt 1^1 salt’ actually needed. In our actual 
implementation, the dimensionless parameters Bz and 
Bi, Eqs. (IT^ and ([5]), are used instead of the chemical 
potentials themselves to simulate the DNA system. The 
values of these parameters for different mixtures of di¬ 
valent and monovalent salts are listed in Sec. IIIII Table 
I. 

Lastly, beside particle addition/deletion moves, one 
also try standard particle translation moves. They are 
carried out exactly like in the case of a canonical Monte- 
Carlo simulation. In a “trial” move from state i to state 
j, an ion is chosen at random and is moved to a random 
position in a volume element surrounding its original po¬ 
sition. The standard Metropolis probability is used for 
the acceptance of such “trial” move: 


/m = min{l, exp[/?(C/i - C/,)]}. (17) 


III. THE SIMULATION MODEL 


fij 

fji 


Bi 

{Niz + l){Ni- + 1) 


exp[/3(C/j 


U,)], 


( 11 ) 


with 


and 


1/2 

Bi = exp(/3^j^ gg^j^) ^3 ^3 , (12) 


1 ^ 1 ,salt = (13) 

is the combined chemical potential of a monovalent salt 
molecule. For a “trial” move where a monovalent salt 
molecule is removed from the system, 

7^ = %^exp[^([/.-t/,)], (14) 

Jji 


Because we are trying to simulate a mixture of salts, to 
improve the system relaxation and to improve the sam¬ 
pling of the system’s phase space, one can also make a 
“trial” move where one Z—ion is added to the system 
and Z monovalent counterions are removed the system. 
For such move, it is easy to show that 


fij 

fji 


BfjV,+ ...(A,+ -Z + l) 

BziNiz + 1) 


exp[/3(t7i - C/j-)]> (15) 


We model the DNA bundle in hexagonal packing as 
a number of DNA molecules arranged in parallel along 
the Z-axis. In the horizontal plane, the DNA molecules 
form a two dimensional hexagonal lattice with lattice 
constant d (the DNA—DNA interaxial distance) (Fig. 
[IJ. Individual DNA molecule is modeled as an impen¬ 
etrable cylinder with negative charges glued on it. The 
charges are positioned in accordance with the locations 
of nucleotide groups along the double-helix structure of 
a B—DNA. The hardcore cylinder has radius of 7A. The 
negative charges are hard spheres of radius 2 A, charge — e 
and lie at a distance of 9A from the DNA axis. This gives 
an averaged DNA radius, r^NA, of Inm. The solvent wa¬ 
ter is treated as a dielectric medium with dielectric con¬ 
stant e = 78 and temperature T = 300° AT. The positions 
of DNA molecules are fixed in space. This mimics the 
constraint on DNA configurational entropy inside viruses 
and other experiments of DNA condensation using diva¬ 
lent counterions in restricted environment. The mobile 
ions in solution are modeled as hard spheres with un¬ 
screened Coulomb interaction (the primitive ion model). 
The colons have radius of cr_ = 2A and charge — e. The 
divalent counterions have radius of az = 2.0, 2.5, or 3.0A 
and charge -|-2e. The interaction between two ions a and 
(3 with radii and charges Qa^p is given by 

jj f QaQp!^‘^aP if ^olP T ^P 

~ \ 00 ii Tap < aa + ap ^ ^ 
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FIG. 1: (Color online) A DNA bundle is modeled as a 
hexagonal lattice with lattice constant d. Individual 
DNA molecule is modeled as a hard-core cylinder with 
negative charges glued on it according to the positions 
of nucleotides of a B—DNA structure. 


Bj 

0.744 X 10® 
2.568 X 10® 
14.48 X 10® 
26.43 X 10® 
56.67 X 10® 
323.82 X 10® 
1302.73 X 10' 


Bl 

0.612 X 10“^ 
0.808 X 10“^ 
1.306 X 10'‘ 
1.580 X 10^ 
2.128 X 10'‘ 
3.715 X lO'^ 
6.601 X 10^ 


cz (mM) 

13.9 ±3.0 

29.9 ±3.4 
74.6 ± 6.2 
99.8 ±5.7 
150.2 ±8.4 
299.6 ± 11.2 
507.1 ± 13.6 


Cl (mM) 
50.0 ±5.6 

50.2 ±4.9 
50.1 ±5.3 

50.3 ±5.4 
50.6 ±6.7 

49.4 ± 6.8 
50.3 ±6.9 


Ph (atm) 
3.183 ±0.001 
4.17 ±0.01 
6.874 ± 0.006 
8.391 ±0.006 
11.42 ±0.02 
20.81 ±0.04 
35.0 ±0.1 


TABLE I: The parameters, and B^, of the salts 
used in the simulation for the reference volume 
~ 2.65 X 10® A® (see text for detail). Columns 3 
and 4 show the corresponding salt concentrations of the 
simulated DNA—free bulk solution. Column 5 shows 
the total pressure of the bulk solutions obtained from 
simulation. 


50 — 2.65 X 10® A®. For a simulation system where 

i is different from 50A, the parameters Bz and Bi are 


/ d 

Bz{d) = B*z ,B,id) = Bl 



( 20 ) 


where = jr^ — r^| is the distance between the ions. 

The simulation is carried out using the periodic bound¬ 
ary condition. Unless explicitly stated, a periodic simu¬ 
lation cell with Nun A = 12 DNA molecules in the hor¬ 
izontal (x, y) plane and 3 full helix periods in the z di¬ 
rection is used. The dimensions of the box are = 3d, 
Ly = 2'/3d and = 102A. This gives, for the volume 
of the simulation box, 

Uceii = 612^3 A^ (19) 

The long-range electrostatic interactions between charges 
in neighboring cells are treated using the Ewald summa¬ 
tion method. In Ref. 00 , it is shown that the macro¬ 
scopic limit is reached when Nona > 7. Our simulation 
cell contains 12 DNA helices, hence it has enough DNA 
molecules to eliminate the finite size effect. Test runs 
with 1, 4, 7 and 12 DNA molecules are carried out to 
verify that this is indeed the case. 

As mentioned above, the DNA bundle is simulated in 
equilibrium with a bulk solution containing two salt con¬ 
centrations: a varying bulk multivalent counterion con¬ 
centrations Cz and a hxed bulk concentration of mono¬ 
valent salt, Cl = 50mM. The detail implementation of 
the GCMC method for this case is described in section 
II. In simulation, the chemical potential of each salt is 
set by fixing the parameters Bi z given by Eq. (Emu). 
In Table HI various values for the parameters B’^ and B* 
that are used in this work for divalent counterion size of 
2A are shown. These values are listed for a reference 
volume that is chosen to have the same dimen¬ 

sions as that of a DNA bundle system with d = 50A, 


In columns 3 and 4 of table I, the resultant salt concen¬ 
trations, Cz and Ci, of the DNA—free solution obtained 
from our GCMC simulations are listed. The divalent salt 
concentration is varied from 14 mM to 507 mM while the 
monovalent salt concentration is kept at approximately 
50 mM. Typical standard deviations in the concentration 
is about 10% in our simulation. This relative error is in 
line with previous GCMC simulations of primitive elec¬ 
trolytes [0. Note that, even though ci is kept constant, 
Bl, (and correspondingly the monovalent salt chemical 
potential galf) ^ constant but actually increases 
with Cz- This is expected because higher cz leads to 
higher free energy cost of adding a monovalent salt to 
the system. 

For each simulation run, about 500-1000 million MC 
moves are carried out depending on the average number 
of ions in the system. To ensure thermalization, about 50 
million initial moves are discarded before doing statistical 
analysis of the result of the simulation. 

In this paper, we are concerned with calculating the 
“effective” DNA—DNA interaction, and correspondingly 
the free energy of assembling DNA bundle. In general, 
this is not a trivial task for a Monte-Carlo simulation 
because the entropy cannot be calculated explicitly. To 
overcome this problem, the Expanded Ensemble method 
0 is implemented. This method allows us to calcu¬ 
late the difference of the system free energies at different 
volumes by sampling these volumes simultaneously in a 
simulation run. By sampling two nearly equal volumes, 
V and U ± AU, and calculate the free energy difference 
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A^l, we can calculate the total pressure of the system: 


P{T,V,M) 


dn{T, V, {/r.}) 


dV 




(21) 

AT/ ^ ^ 


Here are the set of chemical po¬ 

tentials of different ion species. The osmotic pressure of 
the DNA bundle is then obtained by subtracting the total 
pressure of the bulk DNA—free solution, Pb{T, V, {fJ-v}), 
from the total pressure of the DNA system: 


PosmiT, V, {/i,}) = P(T, H, M) - Pb{T, V, {fi,}) (22) 

The total pressure of the bulk solution, Pb{T,V, 
needs to be calculated only once for each set of salt con¬ 
centrations, cz and ci. For reference purpose, their val¬ 
ues are listed in column 5 of Table ID 

All simulations are done using the physics simulation 
library SimEngine develop by one of the author (TTN). 
This library use OpenCL and OpenMP extensions of 
the C programming language to distribute computational 
workloads on multi-core CPU and GPGPU to speed 
up the simulation time. Both molecular dynamics and 
Monte-Carlo simulation methods are supported. In this 
paper the Monte-Carlo module of the library is used. 


IV. RESULT AND DISCUSSION 


A. Counterion mediated DNA DNA interactions 
and the DNA packaging free energy 
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FIG. 2: (Color online) The osmotic pressure of the 
DNA bundle as function of the interaxial DNA distance 
d for different divalent counterion concentration cz 
shown in the inset. The solid lines are guides to the eye. 
The counterion radius is 2.0 A 


volume of the bundle. Per DNA nucleotide base, the 
packaging free energy is given by: 


mdna(^) 


I 




U^'DNA doc 


Posm{d')dV 


I 


N, 


DNA 


( OT T 

Posm{d')^^dd' 
a' 


(23) 


In Fig. [21 the osmotic pressure of DNA bundle at dif¬ 
ferent Cz is plotted as a function of the interaxial DNA 
distance, d for the case the counterion size is 2A. Because 
this osmotic pressure is directly related to the “effective” 
force between DNA molecules at that interaxial distance 
dill, this figure also serves as a plot of DNA—DNA in¬ 
teraction. As one can see, when cz is greater than a value 
around 20mM, there is a short-range attraction between 
two DNA molecules as they approach each other. This 
is the well-known phenomenon of like-charge attraction 
between macroions [ill . MM- It is the result of the 
electrostatic correlations between counterions condensed 
on the surface of each DNA molecule. The attraction 
appears when the distance between these surfaces is of 
the order of the lateral separation between counterions 
(about IdA for divalent counterions). The maximal at¬ 
traction occurs at the distance d ~ 27A, in good agree¬ 
ment with various theoretical and experimental results 
smaller d, the DNA-DNA interaction expe¬ 
riences sharp increase. This can be understood as the 
result of the hardcore repulsion between the counterions. 

From the P-V curve, we can also can calculate the free 
energy, A^p)NA’ packaging DNA into bundles. This 
free energy is nothing but the difference between the free 
energy of a DNA molecule in a bundle and that of an 
individual DNA molecule in the bulk solution (d = oo). 
It can be calculated by integrating the pressure with the 


here I = I.tA is the distance between DNA nucleotides 
along the axis of the DNA. The numerical result for 
dDNA('^*) optimal bundle lattice constant d* is 

plotted in Fig. [3] as function of the cz- Due to the 
limitation of computer simulations, the numerical inte¬ 
gration is performed up to the distance d = 50A only. 
However, this will not change the conclusion of this pa¬ 
per because the omitted integration from d = 50A to 
d = oo only gives an almost constant shift to MdNA ■ 
evident from Fig. |3J the non-monotonic dependence of 
the electrostatic contribution to DNA packaging free en¬ 
ergy is clearly shown. There is an optimal concentration, 
czfi, where the free energy cost of packaging DNA is low¬ 
est. It is negative indicating the tendency of the divalent 
counterions to condense the DNA. At smaller or larger 
concentrations of the counterions, the free energy cost of 
DNA packaging is higher. These results are consistent 
with the correlation theory of DNA reentrant condensa¬ 
tion by multivalent counterions [1, HH and the exper¬ 
iment results on ejecting DNA from bacteriophage under 
varying counterion concentrations Q. However, it must 
be stated, unlike the condensation with counterions of 
higher valence [1, |^, [IJ , the divalent counterions in our 
simulation are not able to decondense the DNA bundle 
within the range of concentration considered. The free 
energy doesnot become positive beyond cz,o- This is in 
line with experimental results 0. 
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FIG. 3: (Color online) The free energy of packaging 
DNA molecules into hexagonal bundles as a function of 
the divalent counterion concentrations. The points are 
results of numerical integration of Posm from Fig. [3l 


Figure [3] gives the short-range attraction among DNA 
molecules to be —O.OdfcsT/base. This is larger than the 
fitted value obtained from the viral DNA ejection exper¬ 
iments M- There are many factors that lead to this 
quantitative discrepancy. Our main approximation is 
that in the simulation, the position of the DNA cylinders 
are straight with infinite bending rigidity. Inside viruses, 
DNA are bent, and the configuration entropy of the DNA 
are not necessary zero, and there is not a perfect hexago¬ 
nal arrangement of DNA cylinder with fixed inter—DNA 
distance. We also neglect the contribution from the re¬ 
gion d > 50A in our integration. The physical parameters 
of the system such as ion sizes, DNA orientations (twist¬ 
ing, frustrations),... [1,[2^|2^ can also affect the strength 
of DNA—DNA short range attraction. All these factors 
are expected to reduce the attraction between the DNA 
compared to our idealized simulation. Nevertheless, the 
non-monotonic electrostatic influence of divalent coun¬ 
terions on DNA-DNA “effective” interaction is clearly 
demonstrated in our idealized simulation. 


B. Role of finite size of counterions 

In all the systems simulated so far, the radius of the 
divalent counterion is fixed at 2.0A. The results agree 
qualitatively and semi-quantitatively with some of the 
experimental results of DNA ejection from capsid with 
MgS 04 salt. However, experimental results also show 
that there is an ion specific effect. There are some sig¬ 
nificant differences in condensations of free DNA, con¬ 
densations of DNA inside viruses when different divalent 
salts such as MgS 04 , MgCl 2 , or MnCl 2 are usedjl, 
This shows that the hydration effect and the entropy of 
the hydrated water molecules are significant and need to 
be properly taken into account when one deals with the 
problem of DNA confinement inside viral capsids. In this 
section, a first step is taken to study this ion specific ef¬ 
fect. Specifically, we study how DNA-DNA interaction 


is affected by changing the radius of the counterions. 

In Fig. m and Fig. [SJ the dependence of DNA- 
DNA ’’effective” interaction on the DNA-DNA separa¬ 
tion distance are plotted for the counterion radii 2.5A 
and SA respectively. Compare to similar plot for the 
case of az = 2.0A (Fig. [2|), we can clearly see that 
the main physics remains when we change the counterion 
size. The DNA-DNA short-range interaction remains ev¬ 
ident. However, the depth and location of the strongest 
attraction change when the counterion size changes. The 
smallest counterions (2A) cause the strongest attraction 
among DNA at smaller distance. This is easily under¬ 
stood, the smaller counterion cause less entropic cost of 
bringing DNA closer to each other. Hence the short- 
range attraction is enhanced. 



FIG. 4: (Color online) The osmotic pressure of the 
DNA bundle as function of the interaxial DNA distance 
d for different divalent counterion concentration cz 
shown in the inset. The solid lines are guides to the eye. 

The counterion radius is 2.5 A 

The change in the equilibrium separation of DNA in 
the bundle is even more evident in Fig. In this fig¬ 
ure, the osmotic pressure (which is proportional to the 
effective DNA-DNA interaction) of the hexagonal DNA 
bundle is plotted as a function of the inter DNA distance 
for three counterion sizes, 2A, 2.5A, and SA, respectively. 
The counterion concentration is chosen to be approxi¬ 
mately I50mM in each simulation. As one can see, the 
first consequence of changing counterion size is obviously 
the equilibrium distance of the DNA bundle. The opti¬ 
mal inter DNA distance, d*, where the short range DNA 
attraction is strongest increases with the counterion ra¬ 
dius. As the counterion radius is increased from 2.0A to 
2 . 5 A to 3 . 0 A d* increases from 26A to 27A then 29A 
respectively. 

However, it is an interesting observation that not only 
the optimum distance d* is shifted by 2az, the interac¬ 
tion between DNA molecules from the distance d* to 00 , 
which is dominated by electrostatics, is shifted by the 
same amount. This is evident as in Fig. [BJd where the 
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FIG. 5: (Color online) The osmotic pressure of the 
DNA bundle as function of the interaxial DNA distance 
d for different divalent counterion concentration cz 
shown in the inset. The solid lines are guides to the eye. 
The counterion radius is az = 3.0A 


horizontal axis for each curve is shifted by 2az- One can 
see that the right side of these curve from the distance 
d* to oo show a good degree of overlapping. 

This is in agreement with the “correlated liquid” na¬ 
ture of DNA—DNA attraction mediated by multivalent 
counterions [mill- In this strongly correlated liquid 
theory of DNA—DNA interaction, the combined system 
of DNA-|-condensed counterions acts as a charged metal¬ 
lic cylinder. The correlations between the condensed 
counterions on the surface of two neighboring DNA in¬ 
duce a short range attraction between them. In this the¬ 
ory, the center of mass of condensed counterion cannot 
approach the DNA surface at a distance less than its ra¬ 
dius, az- Because of this, the effective surface of the 
dressed metallic DNA is lifted off the hare DNA surface 
by a distance of 


X — a z + \ (24) 

where A is the Goy-Chapman length. The length f is half 
the (negative) screening length of the strongly correlated 
liquid of the condensed counterions on the surface of the 
DNA molecule. 


^ 47r(Ze)2 dn 

with fj, the chemical potential of a counterion in the liq¬ 
uid, and n is its two-dimensional density. This screening 
length, 1^1, depends weakly on the ratio, az/roNA- For 
our purpose, it can be considered to be constant. There¬ 
fore, if one considers the correlation-induced attraction 
between two DNA cylinders only works when the clos¬ 
est approach between their surfaces is greater than 2x 
(so that the two DNA’s ’’effective” metallic layers donot 
overlapped), one immediately comes to the conclusion 




FIG. 6: (Color online) a) The osmotic pressure of the 
DNA hexagonal bundle as function of the lattice 
constant, d, for three values of the counterion radius, at 
the same counterion concentration of ISOmM. b) The 
same plot with the horizontal axis shifted by 2az 
showing a good degree of overlapping of the three 
curves with regard to the equilibirum position and the 
attractive electrostatic interaction. 


that the electrostatic like-charged attraction between two 
neighboring DNA cylinders is simply shifted by a dis¬ 
tance of 2az when the radius of the counterion changes. 
This agrees with our simulation results. 

In Fig. [3 the free energy of packaging DNA into an 
hexagonal bundle with the optimal inter—DNA distance, 
d*, is plotted as a function of the counterion concentra¬ 
tions for the three different counterion radii. It can be 
seen clearly that, within the range of counterion concen¬ 
tration studied, there is a quantitative and qualitative 
difference in the free energy of packaging for the three 
sizes of counterion consider. For az = 2A and 2.5A, 
the dependence of the free energy of packaging DNA in 
bundle on the concentration cz is non-monotonic. How¬ 
ever, for the larger counterion size, az = 3A, in the 
range of concentration considered, DNA condense later 
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V. CONCLUSION 



FIG. 7: Free energy per nucleotide base of packaging 
DNA molecules into bundles as a function of counterion 
concentration cz for different radii of the counterions. 


but stronger into hexagonal bundle as the counterion con¬ 
centration increases. This behavior is actually observed 
in experiments. While the non-monotonic behaviors of 
DNA ejection is observed clearly for MgS 04 salt, and 
somewhat evident for MgCl 2 salt, MnCl 2 are known to 
condense DNA in free solution without ever disintegrated 
d, 0 . Our simulation suggests that the difference in the 
hydration radius of the counterions can be used to explain 
such differences. Our results suggests that Mn^+ coun¬ 
terion has larger ion radius. This is in good qualitative 
agreement with computational and EXAFS and X—ray 
studies on divalent counterions hydration shell (see Ta¬ 
ble 3 of reference [13 and the corresponding references 
therein). These work shown that the number of water 
molecules in the hydration shell of ions increases with its 
atomic number. Specifically, as the atomic number of the 
divalent counterions increases from Mg^+, Ca^+, Sr^+ to 
Ba^^, the coordination number increases from 6 to 9 wa¬ 
ter molecules in the hydration shell. Even though, Mn^+ 
hydration was not studied in these works, its atomic num¬ 
ber is higher than that of Ca^+ and Mg^+ ions suggest¬ 
ing that its hydration radius is larger than that of Mg^+ 
counterions. 

It is of importance to note that, according to our 
Fig. [3 although the larger counterions do not produce 
a reentrant non-monotonic behavior, they actually cause 
stronger DNA-DNA attraction energy. Based on what is 
observed from Fig. [3 Fig. 13 Fig. [5] and the horizon¬ 
tally shifted Fig. |B}d, this observation can be explained 
as the result of two effects. First, smaller counterions can 
condense better on DNA, causing a stronger short-range 
like charge attraction among DNA cylinders. However 
they also cause a higher degree of overcharging at larger 
concentrations, so it is costlier to packaging DNA. This 
is evident by the increase in the packaging free energy at 
higher concentration for az = 2A. Secondly, the short- 
range attraction between DNA is shifted to larger d for 
larger counterions. Since one integrates / PdV to find 
the packaging free energy, a simple geometric argument 
shows that the contribution from larger d would dom¬ 
inate this integral, therefore the larger counterions can 
cause lower energy minimum at large concentration. 


In this paper, we use a Grand-Canonical Monte-Carlo 
simulation to study the electrostatics of DNA condensa¬ 
tion, using a primitive model for the screen ions. Specifi¬ 
cally, the effective electrostatic interaction between DNA 
molecules in a hexagonal bundle is computed in the pres¬ 
ence of 50mM monovalent counterions and with varying 
concentration of divalent counterions. The entropy of 
DNA configure fluctuation is suppressed in simulation by 
fixing the position of the DNA cylinders in the bundle. 
Such study can be applied directly to the experimen¬ 
tal problem of DNA ejection from bacteriophages where 
DNA condensed in a strongly confined environment. It is 
shown that, even at the level of non-specific electrostatic 
interaction, divalent counterions can strongly influence 
DNA interaction and packaging. The simulation results 
for divalent counterions with 2.0A radius show that the 
electrostatic free energy of packaging DNA into hexag¬ 
onal bundle varies non-monotonically with the counte¬ 
rion concentration. However, divalent counterions donot 
correlate strong enough with each other to drive DNA 
de-condensation. 

The counterion specificity such as the ion hydration 
radius can influence strongly the qualitative and quan¬ 
titative picture of DNA condensation. Three different 
counterion sizes are studied. They show that the non¬ 
monotonicity changes significantly and disappears as the 
counterion size increases. The most important results of 
this paper are presented in Fig. [Bland Fig. 13 where it is 
shown that increasing counterion radius simply raises the 
” metallic” surface of condensed counterions off the DNA 
and shift the correlation-induced attraction between two 
DNA cylinders by an amount of 2az- This interestingly 
is responsible for making the larger counterions to cause a 
deeper minimum of DNA packaging free energy. In fact, 
in the range of concentration considered in our simulation 
with counterion radius of SA, this free energy keeps go¬ 
ing lower with increasing counterion concentration. Such 
qualitative differences are observed with DNA conden¬ 
sation experiments involving Mg^’*' and Mn^+ counteri¬ 
ons and suggesting that Mn2-|- has bigger ion radius, in 
agreement with previous computational and EXAFS and 
X—ray experimental results. 

Going beyond the scope of DNA ejection experiments, 
we believe the quantitative results of our paper can be 
used to understand many other experiments involving 
DNA and divalent counterions. 
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